#!/bin/bash -e

# based on anno_0216.sh, changed:
# nprot_filter_110416.sh first, instead of last.
function info() {
echo Usage: `basename $0` '-m(somatic) -g(single_cell) -l(interval.bed) snp.vcf indel.vcf'
exit 65
}

while getopts ":w:p:l:mgh" opts
do
case $opts in
	w) win_size=$OPTARG;;
	m) somatic=T;;
	g) single_cell=T;;
	l) interval=$OPTARG;;
	p) out_prefix=$OPTARG;;
	\?) info;;
esac
done
shift $(($OPTIND - 1))

test $# -lt 2 && info


. $var


echo
echo $interval
echo dbsnp version is: $dbsnp_version
echo $snpeff_db_version

# snp cluster remove
if test "$single_cell" = "T"; then
# nprot need snp and indel files
	# nprot_filter_110416.sh $1 $2
    anno_nprot.sh -p$out_prefix $1 $2
	# format_af.sh -p$out_prefix $out_prefix.nprot.filter.vcf
	# anno_info_vcf.pl $out_prefix.nprot.filter.vcf $out_prefix.af.txt > $out_prefix.af.vcf
    
    af_anno.sh -p$out_prefix $1
    
	select_indel_vcf.sh -p$out_prefix.1 $out_prefix.af.vcf && ln -fs $out_prefix.1.indel.vcf $out_prefix.indel.af.vcf
	select_indel_vcf.sh -n -p$out_prefix.1 $out_prefix.af.vcf && ln -fs $out_prefix.1.snp.vcf $out_prefix.snp.af.vcf
else
	
# allele frequency anno
# format_af.sh -p$out_prefix.1 $1
# format_af.sh -p$out_prefix.2 $2 
# anno_info_vcf.pl $1 $out_prefix.1.af.txt > $out_prefix.snp.af.vcf
# anno_info_vcf.pl $2 $out_prefix.2.af.txt > $out_prefix.indel.af.vcf

af_anno.sh -p$out_prefix.snp $1
af_anno.sh -p$out_prefix.indel $2

fi


# snp and indel filter anno
IFS=$'\n'
for dat in $(less $tools_path/script/filters_snp.txt); do

filter=$(echo $dat |cut -f2)
filter_str=$(echo $dat|cut -f1)
filter0=$(echo --filterExpression \"$filter\" --filterName "\"${filter_str}\" ")
filters_snp=$filters_snp$filter0
done
filters_snp=$(echo $filters_snp |perl -pe 's/ $//')

for dat in $(less $tools_path/script/filters_indel.txt); do

filter=$(echo $dat |cut -f2)
filter_str=$(echo $dat|cut -f1)
filter0=$(echo --filterExpression \"$filter\" --filterName "\"${filter_str}\" ")
filters_indel=$filters_indel$filter0
done
filters_indel=$(echo $filters_indel |perl -pe 's/ $//')


echo;echo;echo gatk VariantFiltration
echo java -Xmx$java_memory -jar $gatk -R $ref_genome -T VariantFiltration -o $out_prefix.snp.filtered.vcf --variant $out_prefix.snp.af.vcf $filters_snp |tee .snp.sh

. .snp.sh
	

echo;echo;echo gatk VariantFiltration
echo java -Xmx$java_memory -jar $gatk -R $ref_genome -T VariantFiltration -o $out_prefix.indel.filtered.vcf --variant $out_prefix.indel.af.vcf $filters_indel |tee .indel.sh

. .indel.sh
	
# snp and indel combine
echo;echo;echo gatk CombineVariants
java -Xmx10g $tmp -jar $gatk \
	-R $ref_genome \
	-T CombineVariants \
	--variant $out_prefix.snp.filtered.vcf \
	--variant $out_prefix.indel.filtered.vcf \
	-o $out_prefix.snp.indel.filtered.vcf \
	-genotypeMergeOptions UNSORTED

# snpsift anno
echo;echo;echo snpsift annotate id
java -Xmx10g $tmp -jar $snpsift \
annotate \
-id $data_path/ncbi/dbsnp/All_20150605.vcf.gz \
$out_prefix.snp.indel.filtered.vcf \
> $out_prefix.snpsift.dbsnp.vcf


echo;echo;echo snpsift annotate clinvar
java $tmp -jar $snpsift \
annotate \
-id ${data_path}/ncbi/clinvar/clinvar_20150106.vcf \
$out_prefix.snpsift.dbsnp.vcf \
> $out_prefix.dbsnp.clinvar.vcf


echo;echo;echo snpsift annotate cosmic
java $tmp -jar $snpsift \
annotate \
${data_path}/cosmic/CosmicCodingMuts.vcf \
$out_prefix.dbsnp.clinvar.vcf \
> $out_prefix.id.cosmic.vcf


echo;echo;echo snpsift annotate varType
java $tmp -jar $snpsift \
varType \
$out_prefix.id.cosmic.vcf \
> $out_prefix.cosmic.var_type.snpsift.vcf

:<<var_type
echo;echo;echo gatk VariantAnnotator
java $tmp -jar $gatk \
	-R $ref_genome \
	-T VariantAnnotator \
	--variant $out_prefix.cosmic.var_type.snpsift.vcf \
	-o $out_prefix.cosmic.var_type.vcf \
	-L $interval \
	-A VariantType
var_type

ln -fs $out_prefix.cosmic.var_type.snpsift.vcf $out_prefix.cosmic.var_type.vcf

echo;echo;echo snpsift annotate gwasCat
java $tmp -jar $snpsift \
gwasCat \
-db ${data_path}/snpeff/gwascatalog.txt \
$out_prefix.cosmic.var_type.vcf \
> $out_prefix.var_type.gwas.vcf


echo;echo;echo snpsift dbnsfp
java $tmp -jar $snpsift \
dbnsfp \
-db ${data_path}/dbnsfp/dbNSFP2.9.txt.gz \
-f Ensembl_transcriptid,Uniprot_acc,Interpro_domain,SIFT_score,Polyphen2_HVAR_score \
$out_prefix.var_type.gwas.vcf \
> $out_prefix.gwas.dbnsfp.vcf

# mv $out_prefix.var_type.gwas.vcf $out_prefix.gwas.dbnsfp.vcf
	
# snpeff anno
echo;echo;echo snpeff ann snpeff2gatk
java $tmp -jar $snpeff \
	-c ${snpeff_path}/snpEff.config \
	-v \
	-lof \
	-noStats \
	$snpeff_db_version \
	$out_prefix.gwas.dbnsfp.vcf \
	> $out_prefix.dbnsfp.snpeff2gatk.vcf
	
# vep anno
echo;echo;echo VEP annotation no pick '(can not run in background)'
$tools_path/ensembl-tools-release-82/scripts/variant_effect_predictor/variant_effect_predictor.pl \
--offline \
--force_overwrite \
--dir $data_path/ensembl/ \
--fasta $ref_genome \
--refseq \
--hgvs \
--format vcf \
--vcf \
--input_file $out_prefix.dbnsfp.snpeff2gatk.vcf \
--output_file $out_prefix.vep.vcf


# vcf header modify
echo;echo; echo vcf header modify
cat $out_prefix.vep.vcf \
| awk -F', ' '{if (/##INFO=<ID=\w+, /) { print $$out_prefix","$2","$3","$4}else{print $0}}' \
> $out_prefix.header_modified.vcf

# $bgzip -c $out_prefix.header_modified.vcf > $out_prefix.header_modified.vcf.gz
# $tabix $out_prefix.header_modified.vcf.gz

# high low split
:<<high_low_old
if test "$somatic" = "T"; then
echo filter_high_low.pl
filter_high_low.pl $out_prefix.header_modified.vcf $out_prefix.high.txt $out_prefix.low.txt $out_prefix.header.txt
else
echo filter_pass.pl
filter_pass.pl $out_prefix.header_modified.vcf $out_prefix.high.txt $out_prefix.low.txt $out_prefix.header.txt
fi

cat $out_prefix.header.txt $out_prefix.high.txt > $out_prefix.high.vcf
cat $out_prefix.header.txt $out_prefix.low.txt > $out_prefix.low.vcf

high_low_old

if test test "$somatic" = "T"; then
filter_anno.sh -s -p$out_prefix $out_prefix.header_modified.vcf
else
filter_anno.sh -p$out_prefix $out_prefix.header_modified.vcf
fi

# info extract and format
anno_bt.sh $out_prefix.high.vcf
txt2xls.pl $out_prefix.*multianno.txt $out_prefix.anno.high_confidence.xls
format_jxl.sh $out_prefix.anno.high_confidence.xls $out_prefix.anno.high_confidence.jxl.xls

anno_bt.sh -p$out_prefix.2 $out_prefix.low.vcf
txt2xls.pl $out_prefix.2.*multianno.txt $out_prefix.anno.low_confidence.xls
format_jxl.sh $out_prefix.anno.low_confidence.xls $out_prefix.anno.low_confidence.jxl.xls

:<<single_cell
if test "$single_cell" = "T"; then
	# nprot filter
	nprot_filter.sh $out_prefix.high.vcf

	# impact anno
	impact_anno.sh $out_prefix.nprot.high.vcf
	impact_anno.sh $out_prefix.nprot.low.vcf
else
	impact_anno.sh $out_prefix.high.vcf
fi

single_cell


# impact anno
impact_anno.sh $out_prefix.high.vcf
impact_anno.sh $out_prefix.low.vcf


. $cmd_done

